www.gusucode.com > 精通MATLAB最优化计算全书代码 程序源码 > 随书源码_精通MATLAB最优化计算/第12章 二次规划/ActivdeSet.m

    function [xv,fv] = ActivdeSet(H,c,A,b,x0,AcSet)
sz = size(A);
xx = 1:sz(1);
nonAcSet = zeros(1,1);
m = 1;

for i=1:sz(1)
    if(isempty(find(AcSet == xx(i),1)))
        nonAcSet(m) = i;
        m = m + 1;
    else
        ;
    end
end

invH = inv(H);
bCon = 1;
while bCon
    cl = H*x0 + c;
    Al = A(AcSet,:);
    bl = b(AcSet);
    xm = QuadLagR(H,cl,Al,zeros(length(AcSet),1));

    if xm == 0
        trAl = transpose(Al);
        R = inv(Al*invH*trAl)*Al*invH;
        lamda = R*(H*xm+cl);
        [lmin,index] = min(lamda);
        if lmin>=0
            xv = x0;
            fv = transpose(x0)*H*x0/2+transpose(c)*x0;
            bCon = 0;
            return;
        else
            l = length(AcSet);
            nonAcSet = [nonAcSet AcSet(index)];
            [sa,sb] = sort(nonAcSet);
            nonAcSet = sa;
            tmpAcS = [AcSet(1:(index-1)) AcSet((index+1):l)];
            AcSet = tmpAcS;
        end
    else
        d = xm;
        mAlpha = inf;
        adinAcS = 0;
        for i=1:length(nonAcSet)
            if A(nonAcSet(i),:)*d < 0
                alpha = (b(nonAcSet(i)) - A(nonAcSet(i),:)*x0)/(A(nonAcSet(i),:)*d);
                if alpha < mAlpha
                    mAlpha = alpha;
                    adinAcS = nonAcSet(i);
                end
            end
        end
        mXA = min(1,mAlpha);
        x0 = x0 + mXA*d;
        if mXA < 1
            AcSet = [AcSet adinAcS];
            [sa,sb] = sort(AcSet);
            AcSet = sa;
        else
            cl = H*x0 + c;
            Al = A(AcSet,:);
            bl = b(AcSet);
            trAl = transpose(Al);
            R = inv(Al*invH*trAl)*Al*invH;
            lamda = R*(H*xm+cl);
            [lmin,index] = min(lamda);
            if lmin>=0
                xv = x0;
                fv = transpose(x0)*H*x0/2+transpose(c)*x0;
                bCon = 0;
                return;
            else
                l = length(AcSet);
                nonAcSet = [nonAcSet AcSet(index)];
                [sa,sb] = sort(nonAcSet);
                nonAcSet = sa;
                tmpAcS = [AcSet(1:(index-1)) AcSet((index+1):l)];
                AcSet = tmpAcS;
            end
        end
    end
end